To be able to edit code and run cells, you need to run the notebook yourself. Where would you like to run the notebook?

This notebook takes about 4 minutes to run.

In the cloud (experimental)

Binder is a free, open source service that runs scientific notebooks in the cloud! It will take a while, usually 2-7 minutes to get a session.

On your computer

(Recommended if you want to store your changes.)

  1. Download the notebook:
    notebook.jl
  2. Run Pluto

    (Also see: How to install Julia and Pluto)

  3. Open the notebook file

    Type the saved filename in the open box.

Frontmatter

If you are publishing this notebook on the web, you can set the parameters below to provide HTML metadata. This is useful for search engines and social media.

Author 1

Introduction to Explicit Neutral Models

This notebook illustrates how you can use PMD to optimize network models with an explicit neutral (EN) conductor representation. It will go through the full workflow, consisting of

  • importing OpenDSS network data (and applying transformations as needed);

  • adding OPF-specific data to the model;

  • optimizing the model;

  • inspecting the results.

👀 Reading hidden code
md"""
# Introduction to Explicit Neutral Models

This notebook illustrates how you can use PMD to optimize network models with an explicit neutral (EN) conductor representation. It will go through the full workflow, consisting of
- importing OpenDSS network data (and applying transformations as needed);
- adding OPF-specific data to the model;
- optimizing the model;
- inspecting the results.
"""
19.4 ms

This notebook will make use of the following packages in various places

👀 Reading hidden code
md"""
This notebook will make use of the following packages in various places
"""
245 μs
👀 Reading hidden code
begin
using PowerModelsDistribution
import Ipopt
import DataFrames: DataFrame
end
2.7 s
"/home/runner/.julia/packages/PowerModelsDistribution/Pxoxg/src/.."
pmd_path = joinpath(dirname(pathof(PowerModelsDistribution)), "..")
👀 Reading hidden code
51.4 μs

Building a case

Importing network data

👀 Reading hidden code
13.2 ms

OpenDSS cases with explicit neutral conductors tend to use components in a way which is not supported directly by PMD. These cases model the grounding as a 'reactor', connected between different terminals of the same bus (i.e. from the 4th terminal to ground).

In PMD, a reactor is mapped by default to a line, which will then be connected on both ends to the same bus, which is not allowed. This problem is solved by applying the data model transformation transform_loops!, which will map the reactor to a shunt instead, or merge terminals if the reactor is in fact a short-circuit.

👀 Reading hidden code
355 μs
begin
case_path = pmd_path*"/test/data/en_validation_case_data/test_grounding.dss"
data_eng = parse_file(case_path, transformations=[transform_loops!])
end
👀 Reading hidden code
Circuit has been reset with the 'clear' on line 3 in 'test_grounding.dss'
 
Command 'calcvoltagebases' on line 24 in 'test_grounding.dss' is not supported, skipping.
 
Command 'solve' on line 26 in 'test_grounding.dss' is not supported, skipping.
 
reactors as constant impedance elements is not yet supported, treating reactor.grounding like line
 12.6 s

Note that this test case is very unbalanced. To obtain a more realistic scenario, we reduce the loading by a factor 3.

👀 Reading hidden code
14.6 ms
begin
for (_,load) in data_eng["load"]
load["pd_nom"] *= 1/3
load["qd_nom"] *= 1/3
end
end
👀 Reading hidden code
53.3 μs

Adding OPF specific data

First of all, it is a good practice to remove any bounds that may have been imported from the OpenDSS network data. These can be default bounds, which might not make sense for the specific case and can lead to infeasibility. For example, this occurs for IEEE13, where the default OpenDSS line ratings are too tight for the base case power flows.

👀 Reading hidden code
306 μs
remove_all_bounds!(data_eng)
👀 Reading hidden code
51.0 ms

Voltage bounds

Voltage bounds are naturally expressed between two terminals.

For example, if a wye-connected load is connected between phase a and the neutral n, we want to apply the bounds lb <= |Ua-Un| <= ub, in order to ensure the voltage across the load does not get too low or high. If the neutral is not modeled and assumed to be grounded everywhere, i.e. Un=0, then this can be done equivalently by bounding Ua directtly.

PMD comes with several data model transformations which allow you to quickly apply voltage bounds. The bounds are specified in per unit, relative to the local voltage base (phase-to-ground). If this was not the case, it would be difficult to specify bounds for multiple voltage zones at once.

First consider add_bus_absolute_vbounds!, which allows you to specify absolute bounds, i.e. for each terminal individually. This will set the bus properties vm_lb and vm_ub appropriately.

👀 Reading hidden code
79.2 ms
add_bus_absolute_vbounds!(
data_eng,
phase_lb_pu = 0.8,
phase_ub_pu = 1.2,
neutral_ub_pu = 0.1,
)
👀 Reading hidden code
2.0 s

The ENGINEERING data model includes special properties to easily apply symmetrical bounds for three-phase buses (refer to the documentation for more details). These properties can be populated easily with add_bus_pn_pp_ng_vbounds!, which specifies

  • phase-to-neutral (pn)

  • phase-to-phase (pp)

  • neutral-to-ground (ng)

bounds. Note that since the voltage base is in phase-to-ground, the pp bounds have to multiplied by sqrt(3) for a three-phase network.

👀 Reading hidden code
691 μs
add_bus_pn_pp_ng_vbounds!(
data_eng, [1:3...], 4,
pn_lb_pu = 0.9,
pn_ub_pu = 1.1,
pp_lb_pu = 0.9*sqrt(3),
pp_ub_pu = 1.1*sqrt(3),
ng_ub_pu = 0.1
)
👀 Reading hidden code
24.5 ms

Often, the true goal of voltage bounds is to protect the connected 'units', i.e. loads, generators etc. Therefore, it can suffice to apply voltage bounds only to those buses and terminals which actually have a unit connected to them. This is what add_unit_vbounds! is for. The delta_multiplier is the correction factor which is applied to the specified bounds if a unit has configuration=DELTA.

👀 Reading hidden code
284 μs
add_unit_vbounds!(
data_eng,
lb_pu = 0.91,
ub_pu = 1.09,
delta_multiplier = sqrt(3),
unit_comp_types = ["load"]
)
👀 Reading hidden code
120 ms

Whether it is desirable to apply as many bounds as possible, or only ones which are really needed, will depend on the chosen formulation later on. In any case, these data model transformations make it easy for the user to apply any desired bounds.

You might have noticed that applying all these transformations will lead to redundant constraints. This is resolved in the data model transformation, which determines a set of absolute and pairwise voltage constraints which imply all other ones.

👀 Reading hidden code
400 μs

For example, let's inspect bus b2. This bus has a single-phase load between terminals 1 and 4, so only for that pair the lower bound should be 0.91 instead of 0.9. For an explanation of how the vm_pair_lb property is structured, refer to the documentation.

👀 Reading hidden code
335 μs
begin
data_math_tmp = transform_data_model(data_eng, kron_reduce=false, phase_project=false)
b2_index = data_math_tmp["bus_lookup"]["b2"]
b2_math = data_math_tmp["bus"]["$b2_index"]
b2_math["vm_pair_lb"]
end
👀 Reading hidden code
7.3 s

Adding a generator

So far, we imported a power flow case from OpenDSS. This means that the problem is fully determined, i.e. there is no degrees of freedom left over which to optimize. Therefore, we will add an additional generator to the problem.

👀 Reading hidden code
325 μs
begin
data_eng["generator"] = Dict{String,Any}()
data_eng["generator"]["g1"] = Dict{String,Any}(
"status" => ENABLED,
"bus" => "b2",
"configuration" => WYE,
"connections" => [2,4],
"pg_lb" => [0.0],
"pg_ub" => [20.0],
"qg_lb" => [0.0],
"qg_ub" => [0.0],
"cost_pg_parameters" => fill(0.0, 3)
)
end
👀 Reading hidden code
80.5 μs

Solving a simple OPF problem

👀 Reading hidden code
210 μs

First of all, we apply the data model transformation to go from ENGINEERING to MATHEMATICAL. We need to pass the flags kron_reduce=false and phase_project=false; these activate data model modifications which are required for some formulations, but not for the EN ones.

👀 Reading hidden code
257 μs
data_math = transform_data_model(data_eng, kron_reduce=false, phase_project=false)
👀 Reading hidden code
updated generator 1 cost function with order 3 to a function of order 2: [0.0, 0.0]
 51.1 ms

Before solving the problem, it is important to add initialization values for the voltage variables. Failing to do so will almost always result in solver issues.

For a single-phase equivalent network, this is very easy to do; simply initialize each voltage variable to 1.0+im*0.0, which usually corresponds to the voltage profile without any network load and ignoring all shunts and linecharging (refered to as 'no-load voltage').

The equivalent in a general, multi-conductor network is less trivial; for example, transformers can shift the no-load voltage angle of individual terminals in complicated ways. Therefore, we developed the method add_start_vrvi!, which will infer the no-load voltage for each terminal in the network and add initialization properties to the MATHEMATICAL data model.

👀 Reading hidden code
429 μs
add_start_vrvi!(data_math)
👀 Reading hidden code
830 ms

Now we are ready to solve the OPF problem.

👀 Reading hidden code
244 μs
res = solve_mc_opf(data_math, IVRENPowerModel, Ipopt.Optimizer)
👀 Reading hidden code
❔
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.10, running with linear solver MUMPS 5.5.1.

Number of nonzeros in equality constraint Jacobian...:      428
Number of nonzeros in inequality constraint Jacobian.:      264
Number of nonzeros in Lagrangian Hessian.............:      270

Total number of variables............................:       94
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       22
                     variables with only upper bounds:        0
Total number of equality constraints.................:       76
Total number of inequality constraints...............:       48
        inequality constraints with only lower bounds:       23
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:       25

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.99e-02 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  2.5109835e-04 8.71e-03 5.13e+00  -1.0 7.87e-02  -4.0 1.00e+00 5.64e-01h  1
   2  1.0958071e-04 4.25e-03 2.89e+01  -1.0 1.85e-01  -4.5 1.00e+00 5.12e-01h  1
   3  6.6971358e-05 1.78e-03 4.37e+01  -1.0 3.28e-02  -5.0 1.00e+00 5.81e-01h  1
   4  1.0018333e-04 5.23e-04 4.63e+01  -1.0 2.45e-02  -5.4 1.00e+00 7.07e-01h  1
   5  1.0224534e-04 2.35e-04 2.01e+02  -1.0 4.33e-03  -5.9 1.00e+00 5.50e-01h  1
   6  1.1097993e-04 9.30e-05 4.21e+02  -1.0 3.57e-03  -6.4 1.00e+00 6.05e-01h  1
   7  1.1162713e-04 3.92e-05 1.11e+03  -1.0 7.17e-04  -6.9 1.00e+00 5.78e-01h  1
   8  1.1316823e-04 1.61e-05 2.59e+03  -1.0 6.05e-04  -7.3 1.00e+00 5.89e-01h  1
   9  1.1329233e-04 6.68e-06 6.34e+03  -1.0 1.22e-04  -7.8 1.00e+00 5.85e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.1335869e-04 5.70e-06 3.14e+04  -1.0 1.04e-04  -8.3 1.00e+00 1.47e-01f  3
  11  1.1359999e-04 1.30e-06 1.53e+04  -1.0 5.87e-05  -8.8 1.00e+00 7.71e-01h  1
  12  1.1360176e-04 1.12e-06 1.45e+05  -1.0 1.14e-05  -9.2 1.00e+00 1.35e-01f  3
  13  1.1362502e-04 2.46e-07 7.22e+04  -1.0 5.82e-06  -9.7 1.00e+00 7.74e-01h  1
  14  1.1362759e-04 2.11e-07 6.74e+05  -1.0 4.48e-06 -10.2 1.00e+00 1.38e-01f  3
  15  1.1363741e-04 3.15e-08 2.68e+05  -1.0 2.27e-06 -10.7 1.00e+00 8.12e-01h  1
  16  1.1363742e-04 2.48e-08 2.47e+06  -1.0 4.99e-07 -11.2 1.00e+00 1.62e-01f  3
  17  1.1363823e-04 1.89e-11 4.50e-02  -1.0 1.99e-07 -11.6 1.00e+00 1.00e+00h  1
  18  1.0412787e-04 3.56e-10 4.40e+04  -5.7 3.59e-04 -12.1 9.93e-01 1.00e+00h  1
  19 -2.3291584e-03 6.72e-04 1.69e+04  -5.7 9.19e-02 -12.6 6.80e-01 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20 -3.2606144e-03 2.78e-04 2.17e+03  -5.7 4.48e-02 -13.1 8.71e-01 7.46e-01h  1
  21 -3.2768057e-03 1.97e-04 6.27e+01  -5.7 4.04e-03 -13.5 9.71e-01 7.96e-01h  1
  22 -3.2192508e-03 3.62e-07 5.79e+01  -5.7 2.20e-03 -14.0 7.67e-02 1.00e+00h  1
  23 -3.2230532e-03 4.35e-11 2.82e-04  -5.7 1.38e-04 -14.5 1.00e+00 1.00e+00h  1
  24 -3.2231453e-03 2.59e-14 4.79e-09  -5.7 3.35e-06 -15.0 1.00e+00 1.00e+00h  1
  25 -3.2250651e-03 1.53e-09 6.75e+00  -8.6 7.26e-05 -15.5 1.00e+00 9.63e-01h  1
  26 -3.2250622e-03 9.71e-17 6.41e-11  -8.6 1.02e-07 -15.9 1.00e+00 1.00e+00f  1

Number of Iterations....: 26

                                   (scaled)                 (unscaled)
Objective...............:  -3.2250622204198290e-03   -3.2250622204198290e-03
Dual infeasibility......:   6.4059436658158759e-11    6.4059436658158759e-11
Constraint violation....:   9.7144514654701197e-17    9.7144514654701197e-17
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5108360520094413e-09    2.5108360520094413e-09
Overall NLP error.......:   2.5108360520094413e-09    2.5108360520094413e-09


Number of objective function evaluations             = 35
Number of objective gradient evaluations             = 27
Number of equality constraint evaluations            = 35
Number of inequality constraint evaluations          = 35
Number of equality constraint Jacobian evaluations   = 27
Number of inequality constraint Jacobian evaluations = 27
Number of Lagrangian Hessian evaluations             = 26
Total seconds in IPOPT                               = 1.133

EXIT: Optimal Solution Found.
15.4 s

Inspecting results

👀 Reading hidden code
211 μs

The result dictionary contains the solutions.

👀 Reading hidden code
264 μs
sol_math = res["solution"]
👀 Reading hidden code
12.7 μs

Next, we can transform the solution back to the ENGINEERING data model.

👀 Reading hidden code
231 μs
sol_eng = transform_solution(sol_math, data_math)
👀 Reading hidden code
587 ms

For example, let's inspect the dispatched active generator power.

👀 Reading hidden code
237 μs
sol_eng["generator"]["g1"]["pg"]
👀 Reading hidden code
17.3 μs

The active generator power bounds are not active, 0 <= 5.15 <= 10.0. This means some other bound is active, because since g1 has zero cost, it is normally more optimal to dispatch more active power, which the source bus generator will pay the default price for.

So, let's inspect the voltage at bus b2. We will obtain the voltage base through 'data_math', which will allow us to inspect the voltage in pu. This will make it easier to compare the values against the bounds we specified before.

👀 Reading hidden code
15.3 ms
begin
vbase_b2 = data_math["bus"][string(data_math["bus_lookup"]["b1"])]["vbase"]
v_b2_pu = (sol_eng["bus"]["b2"]["vr"]+im*sol_eng["bus"]["b2"]["vi"])./vbase_b2
vm_b2_pu = abs.(v_b2_pu)
end
👀 Reading hidden code
333 ms

As it turns out, the voltage magnitude constraint om the neutral terminal is binding, vm_b2_pu[4]=0.1.

👀 Reading hidden code
260 μs

Available formulations

👀 Reading hidden code
198 μs

We create a results dictionary which will collect the result dictionary generated by each formulation. We will end this section with a comparison of them.

There are several formulations available which support explicit neutrals. The main distinction is whether the flow variables are current or power variables.

👀 Reading hidden code
277 μs
results = Dict{String,Any}()
👀 Reading hidden code
18.5 μs

Current flow variables

The preferred class of exact formulations are IVR, i.e. with current flow variables (I) and rectangular voltage variables (VR).

IVRENPowerModel is a non-linear formulation.

👀 Reading hidden code
317 μs
results["IVRENPowerModel"] = solve_mc_opf(data_math, IVRENPowerModel, Ipopt.Optimizer)
👀 Reading hidden code
❔
This is Ipopt version 3.14.10, running with linear solver MUMPS 5.5.1.

Number of nonzeros in equality constraint Jacobian...:      428
Number of nonzeros in inequality constraint Jacobian.:      264
Number of nonzeros in Lagrangian Hessian.............:      270

Total number of variables............................:       94
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       22
                     variables with only upper bounds:        0
Total number of equality constraints.................:       76
Total number of inequality constraints...............:       48
        inequality constraints with only lower bounds:       23
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:       25

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.99e-02 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  2.5109835e-04 8.71e-03 5.13e+00  -1.0 7.87e-02  -4.0 1.00e+00 5.64e-01h  1
   2  1.0958071e-04 4.25e-03 2.89e+01  -1.0 1.85e-01  -4.5 1.00e+00 5.12e-01h  1
   3  6.6971358e-05 1.78e-03 4.37e+01  -1.0 3.28e-02  -5.0 1.00e+00 5.81e-01h  1
   4  1.0018333e-04 5.23e-04 4.63e+01  -1.0 2.45e-02  -5.4 1.00e+00 7.07e-01h  1
   5  1.0224534e-04 2.35e-04 2.01e+02  -1.0 4.33e-03  -5.9 1.00e+00 5.50e-01h  1
   6  1.1097993e-04 9.30e-05 4.21e+02  -1.0 3.57e-03  -6.4 1.00e+00 6.05e-01h  1
   7  1.1162713e-04 3.92e-05 1.11e+03  -1.0 7.17e-04  -6.9 1.00e+00 5.78e-01h  1
   8  1.1316823e-04 1.61e-05 2.59e+03  -1.0 6.05e-04  -7.3 1.00e+00 5.89e-01h  1
   9  1.1329233e-04 6.68e-06 6.34e+03  -1.0 1.22e-04  -7.8 1.00e+00 5.85e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.1335869e-04 5.70e-06 3.14e+04  -1.0 1.04e-04  -8.3 1.00e+00 1.47e-01f  3
  11  1.1359999e-04 1.30e-06 1.53e+04  -1.0 5.87e-05  -8.8 1.00e+00 7.71e-01h  1
  12  1.1360176e-04 1.12e-06 1.45e+05  -1.0 1.14e-05  -9.2 1.00e+00 1.35e-01f  3
  13  1.1362502e-04 2.46e-07 7.22e+04  -1.0 5.82e-06  -9.7 1.00e+00 7.74e-01h  1
  14  1.1362759e-04 2.11e-07 6.74e+05  -1.0 4.48e-06 -10.2 1.00e+00 1.38e-01f  3
  15  1.1363741e-04 3.15e-08 2.68e+05  -1.0 2.27e-06 -10.7 1.00e+00 8.12e-01h  1
  16  1.1363742e-04 2.48e-08 2.47e+06  -1.0 4.99e-07 -11.2 1.00e+00 1.62e-01f  3
  17  1.1363823e-04 1.89e-11 4.50e-02  -1.0 1.99e-07 -11.6 1.00e+00 1.00e+00h  1
  18  1.0412787e-04 3.56e-10 4.40e+04  -5.7 3.59e-04 -12.1 9.93e-01 1.00e+00h  1
  19 -2.3291584e-03 6.72e-04 1.69e+04  -5.7 9.19e-02 -12.6 6.80e-01 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20 -3.2606144e-03 2.78e-04 2.17e+03  -5.7 4.48e-02 -13.1 8.71e-01 7.46e-01h  1
  21 -3.2768057e-03 1.97e-04 6.27e+01  -5.7 4.04e-03 -13.5 9.71e-01 7.96e-01h  1
  22 -3.2192508e-03 3.62e-07 5.79e+01  -5.7 2.20e-03 -14.0 7.67e-02 1.00e+00h  1
  23 -3.2230532e-03 4.35e-11 2.82e-04  -5.7 1.38e-04 -14.5 1.00e+00 1.00e+00h  1
  24 -3.2231453e-03 2.59e-14 4.79e-09  -5.7 3.35e-06 -15.0 1.00e+00 1.00e+00h  1
  25 -3.2250651e-03 1.53e-09 6.75e+00  -8.6 7.26e-05 -15.5 1.00e+00 9.63e-01h  1
  26 -3.2250622e-03 9.71e-17 6.41e-11  -8.6 1.02e-07 -15.9 1.00e+00 1.00e+00f  1

Number of Iterations....: 26

                                   (scaled)                 (unscaled)
Objective...............:  -3.2250622204198290e-03   -3.2250622204198290e-03
Dual infeasibility......:   6.4059436658158759e-11    6.4059436658158759e-11
Constraint violation....:   9.7144514654701197e-17    9.7144514654701197e-17
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5108360520094413e-09    2.5108360520094413e-09
Overall NLP error.......:   2.5108360520094413e-09    2.5108360520094413e-09


Number of objective function evaluations             = 35
Number of objective gradient evaluations             = 27
Number of equality constraint evaluations            = 35
Number of inequality constraint evaluations          = 35
Number of equality constraint Jacobian evaluations   = 27
Number of inequality constraint Jacobian evaluations = 27
Number of Lagrangian Hessian evaluations             = 26
Total seconds in IPOPT                               = 0.027

EXIT: Optimal Solution Found.
32.3 ms

IVRQuadraticENPowerModel is an equivalent quadratic formulation.

👀 Reading hidden code
224 μs
results["IVRQuadraticENPowerModel"] = solve_mc_opf(data_math, IVRQuadraticENPowerModel, Ipopt.Optimizer)
👀 Reading hidden code
❔
This is Ipopt version 3.14.10, running with linear solver MUMPS 5.5.1.

Number of nonzeros in equality constraint Jacobian...:      335
Number of nonzeros in inequality constraint Jacobian.:      232
Number of nonzeros in Lagrangian Hessian.............:      212

Total number of variables............................:      103
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       23
                     variables with only upper bounds:        0
Total number of equality constraints.................:       86
Total number of inequality constraints...............:       44
        inequality constraints with only lower bounds:       21
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:       23

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.99e-02 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -3.5820658e-03 4.38e-03 1.50e+01  -1.0 2.55e-01  -4.0 5.58e-01 7.81e-01H  1
   2 -4.1531552e-04 1.10e-03 3.67e+01  -1.0 1.13e-01  -4.5 1.97e-01 9.88e-01h  1
   3 -9.7936762e-05 2.38e-05 4.51e+01  -1.0 1.27e-02  -5.0 1.00e+00 9.90e-01h  1
   4  1.0952651e-04 4.87e-06 1.48e-01  -1.0 7.85e-03  -5.4 1.00e+00 1.00e+00h  1
   5  1.0064762e-04 9.21e-09 3.02e+03  -5.7 3.27e-04  -5.9 9.93e-01 1.00e+00h  1
   6 -2.2532511e-03 6.29e-04 9.43e+02  -5.7 8.89e-02  -6.4 6.87e-01 1.00e+00f  1
   7 -3.2607442e-03 2.78e-04 1.26e+02  -5.7 4.81e-02  -6.9 8.67e-01 7.56e-01h  1
   8 -3.2765187e-03 1.96e-04 3.46e+00  -5.7 4.02e-03  -7.3 9.73e-01 7.97e-01h  1
   9 -3.2201531e-03 3.58e-07 3.17e+00  -5.7 2.15e-03  -7.8 8.33e-02 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -3.2230792e-03 1.03e-09 1.83e-04  -5.7 1.06e-04  -8.3 1.00e+00 1.00e+00h  1
  11 -3.2231450e-03 5.17e-13 2.10e-09  -5.7 2.40e-06  -8.8 1.00e+00 1.00e+00h  1
  12 -3.2250651e-03 1.53e-09 6.07e+00  -8.6 7.27e-05  -9.2 9.96e-01 9.63e-01h  1
  13 -3.2250621e-03 1.00e-15 2.39e-10  -8.6 1.07e-07  -9.7 1.00e+00 1.00e+00f  1

Number of Iterations....: 13

                                   (scaled)                 (unscaled)
Objective...............:  -3.2250621414144703e-03   -3.2250621414144703e-03
Dual infeasibility......:   2.3880537303926607e-10    2.3880537303926607e-10
Constraint violation....:   9.9974282324888364e-16    9.9974282324888364e-16
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5086454645590676e-09    2.5086454645590676e-09
Overall NLP error.......:   2.5086454645590676e-09    2.5086454645590676e-09


Number of objective function evaluations             = 15
Number of objective gradient evaluations             = 14
Number of equality constraint evaluations            = 15
Number of inequality constraint evaluations          = 15
Number of equality constraint Jacobian evaluations   = 14
Number of inequality constraint Jacobian evaluations = 14
Number of Lagrangian Hessian evaluations             = 13
Total seconds in IPOPT                               = 0.015

EXIT: Optimal Solution Found.
6.7 s

Reduced models only create explicit series current variables, and create the total current variables as linear expressions of those. Since branches tend to be the dominate component in number, this can lead to a big reduction in the number of variables.

IVRReducedENPowerModel is the branch-reduced version of IVRENPowerModel.

👀 Reading hidden code
365 μs
results["IVRReducedENPowerModel"] = solve_mc_opf(data_math, IVRReducedENPowerModel, Ipopt.Optimizer)
👀 Reading hidden code
❔
This is Ipopt version 3.14.10, running with linear solver MUMPS 5.5.1.

Number of nonzeros in equality constraint Jacobian...:      364
Number of nonzeros in inequality constraint Jacobian.:      264
Number of nonzeros in Lagrangian Hessian.............:      270

Total number of variables............................:       62
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       22
                     variables with only upper bounds:        0
Total number of equality constraints.................:       44
Total number of inequality constraints...............:       48
        inequality constraints with only lower bounds:       23
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:       25

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.99e-02 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  2.5109816e-04 8.71e-03 5.13e+00  -1.0 7.87e-02  -4.0 1.00e+00 5.64e-01h  1
   2  1.0958060e-04 4.25e-03 2.89e+01  -1.0 1.85e-01  -4.5 1.00e+00 5.12e-01h  1
   3  6.6971301e-05 1.78e-03 4.37e+01  -1.0 3.28e-02  -5.0 1.00e+00 5.81e-01h  1
   4  1.0018331e-04 5.23e-04 4.63e+01  -1.0 2.45e-02  -5.4 1.00e+00 7.07e-01h  1
   5  1.0224533e-04 2.35e-04 2.01e+02  -1.0 4.33e-03  -5.9 1.00e+00 5.50e-01h  1
   6  1.1097993e-04 9.30e-05 4.21e+02  -1.0 3.57e-03  -6.4 1.00e+00 6.05e-01h  1
   7  1.1162713e-04 3.92e-05 1.11e+03  -1.0 7.17e-04  -6.9 1.00e+00 5.78e-01h  1
   8  1.1316823e-04 1.61e-05 2.59e+03  -1.0 6.05e-04  -7.3 1.00e+00 5.89e-01h  1
   9  1.1329233e-04 6.68e-06 6.34e+03  -1.0 1.22e-04  -7.8 1.00e+00 5.85e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.1335869e-04 5.70e-06 3.14e+04  -1.0 1.04e-04  -8.3 1.00e+00 1.47e-01f  3
  11  1.1359999e-04 1.30e-06 1.53e+04  -1.0 5.87e-05  -8.8 1.00e+00 7.71e-01h  1
  12  1.1360176e-04 1.12e-06 1.45e+05  -1.0 1.14e-05  -9.2 1.00e+00 1.35e-01f  3
  13  1.1362502e-04 2.46e-07 7.22e+04  -1.0 5.82e-06  -9.7 1.00e+00 7.74e-01h  1
  14  1.1362759e-04 2.11e-07 6.74e+05  -1.0 4.48e-06 -10.2 1.00e+00 1.38e-01f  3
  15  1.1363741e-04 3.15e-08 2.68e+05  -1.0 2.27e-06 -10.7 1.00e+00 8.12e-01h  1
  16  1.1363742e-04 2.48e-08 2.47e+06  -1.0 4.99e-07 -11.2 1.00e+00 1.62e-01f  3
  17  1.1363823e-04 9.70e-12 4.50e-02  -1.0 1.99e-07 -11.6 1.00e+00 1.00e+00h  1
  18  1.0412787e-04 3.56e-10 4.40e+04  -5.7 3.59e-04 -12.1 9.93e-01 1.00e+00h  1
  19 -2.3291583e-03 6.72e-04 1.41e+04  -5.7 9.19e-02 -12.6 6.80e-01 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20 -3.2606144e-03 2.78e-04 1.81e+03  -5.7 4.48e-02 -13.1 8.71e-01 7.46e-01h  1
  21 -3.2768057e-03 1.97e-04 5.24e+01  -5.7 4.04e-03 -13.5 9.71e-01 7.96e-01h  1
  22 -3.2192508e-03 3.62e-07 4.84e+01  -5.7 2.20e-03 -14.0 7.67e-02 1.00e+00h  1
  23 -3.2230532e-03 4.35e-11 2.82e-04  -5.7 1.38e-04 -14.5 1.00e+00 1.00e+00h  1
  24 -3.2231453e-03 2.59e-14 4.79e-09  -5.7 3.35e-06 -15.0 1.00e+00 1.00e+00h  1
  25 -3.2250651e-03 1.53e-09 6.76e+00  -8.6 7.26e-05 -15.5 1.00e+00 9.63e-01h  1
  26 -3.2250622e-03 8.67e-17 6.41e-11  -8.6 1.02e-07 -15.9 1.00e+00 1.00e+00f  1

Number of Iterations....: 26

                                   (scaled)                 (unscaled)
Objective...............:  -3.2250622204214210e-03   -3.2250622204214210e-03
Dual infeasibility......:   6.4059172000596807e-11    6.4059172000596807e-11
Constraint violation....:   8.6736173798840355e-17    8.6736173798840355e-17
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5108425993955030e-09    2.5108425993955030e-09
Overall NLP error.......:   2.5108425993955030e-09    2.5108425993955030e-09


Number of objective function evaluations             = 35
Number of objective gradient evaluations             = 27
Number of equality constraint evaluations            = 35
Number of inequality constraint evaluations          = 35
Number of equality constraint Jacobian evaluations   = 27
Number of inequality constraint Jacobian evaluations = 27
Number of Lagrangian Hessian evaluations             = 26
Total seconds in IPOPT                               = 0.022

EXIT: Optimal Solution Found.
4.6 s

IVRReducedQuadraticENPowerModel is the branch-reduced version of IVRQuadraticENPowerModel.

👀 Reading hidden code
269 μs
results["IVRReducedQuadraticENPowerModel"] = solve_mc_opf(data_math, IVRReducedQuadraticENPowerModel, Ipopt.Optimizer)
👀 Reading hidden code
❔
This is Ipopt version 3.14.10, running with linear solver MUMPS 5.5.1.

Number of nonzeros in equality constraint Jacobian...:      271
Number of nonzeros in inequality constraint Jacobian.:      232
Number of nonzeros in Lagrangian Hessian.............:      212

Total number of variables............................:       71
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       23
                     variables with only upper bounds:        0
Total number of equality constraints.................:       54
Total number of inequality constraints...............:       44
        inequality constraints with only lower bounds:       21
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:       23

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.99e-02 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -3.5820661e-03 4.38e-03 1.50e+01  -1.0 2.55e-01  -4.0 5.58e-01 7.81e-01H  1
   2 -4.1531578e-04 1.10e-03 3.67e+01  -1.0 1.13e-01  -4.5 1.97e-01 9.88e-01h  1
   3 -9.7936858e-05 2.38e-05 4.51e+01  -1.0 1.27e-02  -5.0 1.00e+00 9.90e-01h  1
   4  1.0952651e-04 4.87e-06 1.48e-01  -1.0 7.85e-03  -5.4 1.00e+00 1.00e+00h  1
   5  1.0064762e-04 9.21e-09 3.02e+03  -5.7 3.27e-04  -5.9 9.93e-01 1.00e+00h  1
   6 -2.2532513e-03 6.29e-04 9.43e+02  -5.7 8.89e-02  -6.4 6.87e-01 1.00e+00f  1
   7 -3.2607442e-03 2.78e-04 1.26e+02  -5.7 4.81e-02  -6.9 8.67e-01 7.56e-01h  1
   8 -3.2765187e-03 1.96e-04 3.46e+00  -5.7 4.02e-03  -7.3 9.73e-01 7.97e-01h  1
   9 -3.2201531e-03 3.58e-07 3.17e+00  -5.7 2.15e-03  -7.8 8.33e-02 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -3.2230792e-03 1.03e-09 1.83e-04  -5.7 1.06e-04  -8.3 1.00e+00 1.00e+00h  1
  11 -3.2231450e-03 5.17e-13 2.10e-09  -5.7 2.40e-06  -8.8 1.00e+00 1.00e+00h  1
  12 -3.2250651e-03 1.53e-09 6.75e+00  -8.6 7.27e-05  -9.2 1.00e+00 9.63e-01h  1
  13 -3.2250621e-03 9.99e-16 6.28e-11  -8.6 1.07e-07  -9.7 1.00e+00 1.00e+00f  1

Number of Iterations....: 13

                                   (scaled)                 (unscaled)
Objective...............:  -3.2250621417804770e-03   -3.2250621417804770e-03
Dual infeasibility......:   6.2804160871431211e-11    6.2804160871431211e-11
Constraint violation....:   9.9898388172814379e-16    9.9898388172814379e-16
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5082794642103833e-09    2.5082794642103833e-09
Overall NLP error.......:   2.5082794642103833e-09    2.5082794642103833e-09


Number of objective function evaluations             = 15
Number of objective gradient evaluations             = 14
Number of equality constraint evaluations            = 15
Number of inequality constraint evaluations          = 15
Number of equality constraint Jacobian evaluations   = 14
Number of inequality constraint Jacobian evaluations = 14
Number of Lagrangian Hessian evaluations             = 13
Total seconds in IPOPT                               = 0.012

EXIT: Optimal Solution Found.
4.6 s

Power flow variables

We also included a single formulation with power flow variables. Because it also has rectangular voltage variables and legacy reasons, it is referred to as ACRENPowerModel.

👀 Reading hidden code
318 μs
results["ACRENPowerModel"] = solve_mc_opf(data_math, ACRENPowerModel, Ipopt.Optimizer)
👀 Reading hidden code
❔
This is Ipopt version 3.14.10, running with linear solver MUMPS 5.5.1.

Number of nonzeros in equality constraint Jacobian...:     1672
Number of nonzeros in inequality constraint Jacobian.:      232
Number of nonzeros in Lagrangian Hessian.............:     1650

Total number of variables............................:       61
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       23
                     variables with only upper bounds:        0
Total number of equality constraints.................:       60
Total number of inequality constraints...............:       44
        inequality constraints with only lower bounds:       21
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:       23

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.60e+06 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -8.7121238e+02 3.17e+04 5.77e+16  -1.0 1.60e+04    -  6.75e-01 9.80e-01f  1
   2 -9.3161445e+00 3.27e+02 1.45e+16  -1.0 1.57e+04    -  8.71e-01 9.90e-01h  1
   3 -7.9560572e-02 2.63e+00 3.20e+15  -1.0 1.68e+02    -  1.00e+00 9.92e-01h  1
   4 -3.1993878e-03 9.89e-05 4.91e+14  -1.0 1.35e+00    -  1.00e+00 1.00e+00h  1
   5 -3.7163678e-03 1.57e-05 5.52e+13  -1.0 1.04e-02    -  1.00e+00 1.00e+00h  1
   6 -4.2146023e-03 1.45e-05 6.47e+12  -1.0 1.02e-02    -  1.00e+00 1.00e+00h  1
   7 -4.2406677e-03 4.15e-08 2.88e+10  -1.0 4.98e-04    -  1.00e+00 1.00e+00h  1
   8 -4.2406237e-03 2.42e-08 2.15e+05  -1.0 1.06e-06    -  1.00e+00 1.00e+00h  1
   9 -4.2406189e-03 2.37e-08 4.63e-05  -1.0 3.91e-08    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -4.4181425e-03 1.85e-06 3.06e+05  -8.6 3.61e-03    -  9.69e-01 1.00e+00f  1
  11 -1.4242374e-02 5.66e-03 1.02e+05  -8.6 1.99e-01    -  6.65e-01 1.00e+00f  1
  12 -1.7214514e-02 3.29e-03 1.47e+04  -8.6 8.31e-02    -  8.56e-01 5.41e-01h  1
  13 -1.6993129e-02 2.44e-05 1.57e+03  -8.6 1.18e-02    -  8.93e-01 1.00e+00h  1
  14 -1.6991466e-02 4.51e-08 6.40e-01  -8.6 8.97e-05    -  1.00e+00 1.00e+00h  1
  15 -1.6991475e-02 1.86e-08 7.53e-09  -8.6 4.21e-08    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 15

                                   (scaled)                 (unscaled)
Objective...............:  -1.6991475049121411e-02   -1.6991475049121411e-02
Dual infeasibility......:   7.5255744949028790e-09    7.5255744949028790e-09
Constraint violation....:   1.2076003556082248e-14    1.8626451492309570e-08
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5059561617557849e-09    2.5059561617557849e-09
Overall NLP error.......:   2.5059561617557849e-09    1.8626451492309570e-08


Number of objective function evaluations             = 16
Number of objective gradient evaluations             = 16
Number of equality constraint evaluations            = 16
Number of inequality constraint evaluations          = 16
Number of equality constraint Jacobian evaluations   = 16
Number of inequality constraint Jacobian evaluations = 16
Number of Lagrangian Hessian evaluations             = 15
Total seconds in IPOPT                               = 1.019

EXIT: Optimal Solution Found.
9.6 s

We mentioned before that the IVR formulations are preferred for EN models. This is because in EN models, the voltage magnitude cannot be bounded below for some terminals (the ones belonging to the neutral conductor). In a formulation with power flow variables, this means that KCL cannot be enforced for those terminals.

In short, ACR allows non-physical groundings of the neutral conductor, and is therefore a relaxation of the original problem. To demonstrate this, let's create a summary of the obtained objective values and generator set points.

👀 Reading hidden code
336 μs
formulationobjective valueg1 pg
StringFloat64Float64
1
"ACRENPowerModel"
-0.0169915
20.0
2
"IVRENPowerModel"
-0.00322506
5.07975
3
"IVRQuadraticENPowerModel"
-0.00322506
5.07975
4
"IVRReducedENPowerModel"
-0.00322506
5.07975
5
"IVRReducedQuadraticENPowerModel"
-0.00322506
5.07975
begin
forms = sort([keys(results)...])
sol_engs = Dict(f=>transform_solution(results[f]["solution"], data_math) for f in forms)
DataFrame(
"formulation" => forms,
"objective value" => [results[f]["objective"] for f in forms],
"g1 pg" => [sol_engs[f]["generator"]["g1"]["pg"][1] for f in forms],
)
end
👀 Reading hidden code
138 ms

This table illustrates that ACRENPowerModel is a relaxation of the other ones. However, note that it is not guaranteed that the objective value will be lower, because all problems are only solved to local optimality. And in fact, the ACR solution is very sensitive to changes in the initialization. The optional virtual groundings seem to introduce many potential local optima.

👀 Reading hidden code
284 μs